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Abstract In this work we present a very simple and efficient numerical scheme 
which can be applied to study the dynamics of bosonic systems like, for in¬ 
stance, spinor Bose-Einstein condensates with nonlocal interactions but equally 
well works for Fermi gases. The method we use is a modification of well known 
Split Operator Method (SOM). We carefully examine this algorithm in the 
case of F = 1 spinor Bose-Einstein condensate without and with dipolar inter¬ 
actions and for strongly interacting two-component Fermi gas. Our extension 
of the SOM method has many advantages: it is fast, stable, and keeps constant 
all the physical constraints (constants of motion) at high level. 

Keywords Split Operator Method ■ SOM • Gross-Pitaevskii equation ■ GP • 
spinor condensate • dipolar condensate • Bose-Einstein condensate • BEC • 
degenerate Fermi gas • nonlinear partial integro-differential set of equations • 
PDE • NoPDE • PIDE 


1 Introduction 

The first dilute atomic Bose-Einstein condensates (BEC) were realized exper¬ 
imentally in 1995 by the groups of: E. Cornell and C. Wieman for rubidium 
PQ; W. Ketterle for sodium [2j, and R. Hulet for lithium atoms [3]. In these 
experiments, magnetic moments (associated with the electrons’ and nucleus’ 
spins) of very low energy atoms followed the external trapping magnetic field. 
Because of spin-spin interactions the projections of magnetic moments can 
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change and thus atoms are no longer kept by the magnetic trap and they es¬ 
cape, so eventually only one component gas remains trapped (with frozen spin 
degree of freedom), and is described by the scalar wave function [4] . 

After experiments with optical traps the atoms’ spin degree of freedom is 
not constrained to a single component only [5j what allows to study a spin 
dynamics due to interactions. But in this case the condensate wave function 
is no longer a single scalar function. It has now 2 F + 1 components describing 
condensates of atoms of all possible spin projections of the total atom spin F. 
Such a system is known as a spinor condensate [6]. 

Soon after the Bose-Einstein condensate was achieved experimentalists suc¬ 
cessfully cooled atomic Fermi gas below the degeneracy temperature [7]. Since 
at very low temperatures the contact interactions in a single-component spin- 
polarized Fermi gas are excluded by the statistics, the experimental realization 
of quantum degeneracy requires trapping of two kinds of atoms. Therefore, in 
the case of fermions the system also becomes multicomponent. In this work 
we describe an efficient method for solving the dynamics of both spinor con¬ 
densates and two-component Fermi gases. 

In addition to the short-range interactions between atoms we include in our 
analysis also the long-range dipolar interactions between magnetic moments of 
atoms. To study the dynamics of such a system of cold bosonic atoms we ap¬ 
ply the mean-field approximation. For fermions, a hydrodynamical approach 
including the gradient corrections [H] to the standard Thomas-Fermi approxi¬ 
mation 0, followed by the inverse Madelung transformation m is employed. 
From the numerical point of view it means that we must solve the system 
of nonlinear partial integro-differential set of equations. This is a demanding 
task. Fortunately, an advantage is that the set of equations under consider¬ 
ation has some constants of motion. Then the algorithm which is developed 
can be verified against how well these quantities are preserved during the evo¬ 
lution. Of course, there exist algorithms belonging to the very fast developing 
group of structure-preserving numerical methods (called often as geometric 
numerical integration methods) which conserve such quantities inherently, by 
construction mm- However, these algorithms, to the best of our knowledge, 
were developed just for ordinary differential equations nn or for only certain 
partial differential equations mainly in a one-dimensional space ini. Our task 
is obviously much more complex. In such a case the procedure of checking the 
conservation of constants of motion is the only one (the other could be the 
comparison with existing analytical solutions which is not, however, our case) 
enhancing our confidence in the correctness of the algorithm. 

The paper is organized as follows. We start with bosons in Section^ where 
we introduce the Gross-Pitaevskii (GP) equation for the scalar condensate 
and extend it to a spinor BEC. In Section [3] we briefly describe a well known 
version of SOM for the case of a single GP equation. Next, we generalize 
this method to the spinor case (Sections U) and study its performance, in 
particular the conservation of constants of motion without dipolar (Section 
[5]) and including dipolar interactions (Sections 16171) . Then we discuss the way 
we treat the system of degenerate fermionic atoms (Section [5]). We show that 
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within the widely accepted approximations for Fermi systems their dynamics 
can be efficiently studied with the help of SOM method. We conclude in section 

M 


2 Description of a condensate — from the second quantization to 
the mean-field approximation 


We start with bosonic systems. The Heisenberg equation of motion written in 
the second quantization formalism is 

= (^o + H c + Hd) $ , (1) 

ot 

where ^(r) is a bosonic field operator, which annihilates an atom at position r. 
The single-particle Hamiltonian Hq describes the contribution from the kinetic 
and potential energies and equals 

#o = -|-V 2 + V trap , (2) 

2 m 

where m is a mass of an atom, Vt ra p = |w 2 (i 2 + y 2 + z 2 ) is a trapping po¬ 
tential (without any loss of generality we assume here a spherically-symmetric 
trap). H c is the Hamiltonian related to the contact interactions which for a 
scalar condensate equals 

H c = g$(r)'4>(r) ( 3 ) 

with the constant g = Ai:h 2 a s /m characterizing the atom-atom interaction (a s 
is the scattering length - the parameter which is sufficient to describe the low 
temperature collisions HD- 

In the case of a condensate of F = 1 atoms one has the spinor field operator 

V’(r) = (^i( r ),V' 0 ( r ),'(/’-i( r )) T , (4) 

where the field operator ^ (r) annihilates an atom in the hyperfine state IF, i) 
at point r. Now, the interaction part of the Hamiltonian in Eq. m readqj 

H c = coi’l (r)V-i(r) + c 2 (r)Fy^-(r)) F, (5) 


where Co = 47rTi 2 (2a2 + ao)/3m and c 2 = 47r/i 2 (a 2 — ao)/3m determine the 
strength of spin-preserving and spin-changing collisions, respectively Q~3] (with 
ao ,2 being the s-wave scattering lengths for the total spin of colliding atoms 
equal to 0 and 2, respectively 0), and F is a spin-one vector built of F = 1 
spin matrices. We use standard definition of F = (F x , F y ,F z ) with 


F x 


h 

71 



J7 — 
> r y ~ 


h 

71 


/0-i 0\ (1 0 0 \ 

i 0 -i , F z = h 0 0 0 

\0 i 0 / \0 0 -1 J 


( 6 ) 


1 


Summation over repeated indices is assumed. 
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The term H c is responsible for the spin dynamics, using matrix notation it 
equals 

(H cll H c 10 0 \ 

H c = H* 10 H c00 H c0 _ i , (7) 

V 0 J 

with definitions 

H c n = (co + c 2 )^l^i + (c 0 + c 2 )ipl'ipo + (co - c 2 )V’li^-i (8a) 

Hcoo = (co + + cq^I^o + (co + c 2 )V>l 1 'f/’-i (8b) 

H c - 1-1 = (c 0 - c 2 )^|^i + (c 0 + C2)i>li >0 + (c 0 + c 2 )^l 1 ^_i (8c) 


Hcio = c 2 ^Li^o , H c 0 _i = c 2 ^J'0i 


= 0 . 


( 8 d) 


Finally, the third term in Eq. m is related to the dipolar interactions and 
is written as 

H d = J d 3 r'ft(r')V d {r-r')$(r'), (9) 

where the interaction energy of two atoms with magnetic moments 7 F 1 and 
7 F 2 (7 is a gyromagnetic coefficient), positioned at r and r' is 


Vd{ r,r') = 7 s 


F 1 F 2 


-3r 


[F, (r-r')] [F a (r - r')] 


( 10 ) 


As can be seen from GOD the spin projection of a pair of colliding atoms 
can change at most by 2 , while the spin projection of a single atom changes 
maximally by 1. Therefore, the matrix H d gets tridiagonal with nonvanishing 
elements on diagonal equal to H daa = aH d n and off-diagonal elements given 
by H dat<x -1 = vt 4 - «)(3 + a)/12 H dW (a = 1,0, -1). 

The equation of motion in a matrix form looks as follows 



H d n 


dlO 


Hcio 

Fod 


H, 


c0—1 


r H d 10 
He 00 

+ Hi 10 H 0 


0 

H c o-i + H d 10 

+ H c - i_i — H d n 



( 11 ) 


Unfortunately, we can not solve the above set of equations because of their 
operator character. One needs additional approximations. From now on we will 
totally neglect the quantum fluctuations and make the following substitution 
for the held operators V>(r) (and V’l(r)) 

VKr) ->■ V’M, 


(12) 
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where r) is a complex function which we often call an order parameter pil l 11 1 . 
Obviously, this is only an approximation^ Nevertheless this kind of treatment 
is justified |15] and agrees well with many experiments [16] . In the case of a 
spinor condensate we have 

(V'i(r),V>o(r),V>-i(r)) T -t (t/>i(r),t/) 0 (r),t/>_i(r)) T , (13) 


where ipi( r), i = 1,0, — 1 are complex wave functions. 

The Heisenberg equation © with an additional condition (fT21) leads us to 
the well known Gross-Pitaevskii equation E supplemented by the nonlocal 
term due to the dipolar interactions 


+ Vtrap(r) + ffIV’WI 2 


+ n 2 7 2 


dV 


r'|3 


-3 


{z-z'f 


r'|5 


IV'(r ')| 2 V'W- 


(14) 


It is worth to notice that the above equation has the form of the Schrodinger 
equation with additional nonlinear and nonlocal terms which take into ac¬ 
count the mean-field and dipolar energies of interacting bosons. For the spinor 
condensates one can obtain 


0 / M / M 

lh dt ( j = (#0 + H c + H d ) I f/’o j • (15) 


One component Gross-Pitaevskii equation without and with additional dipo¬ 
lar interactions, Eq. m, has been already attempted by several authors [18]. 
Also, the spinor version of the GP equation, Eq. (usd, but without including 
dipole-dipole term H d has been studied m- Here, for the first time we in¬ 
vestigate the quality of a numerical algorithm while atomic transfer between 
different spin components, due to dipolar interactions, is allowed. Introducing 
the dipolar interactions into the spinor condensate is by no means a trivial 
extention. In fact, the resulting atomic flow between components corresponds 
to the famous Einstein-de Haas effect [2D]. From the numerical point of view 
it means that another constant of motion emerges in this case, which is the 
projection of the total angular momentum. Any proposed algorithm should be 
checked for the quality of the conservation of this new constant of motion. 

In the rest of this paper we will be often using the oscillatory units, in 
which the units of time, length, energy, coupling constant, and gyromagnetic 

coefficient are given by r = 1/w, a = yj^ , E = Hu, Ea 3 , and y/woP/h, 
respectively. 

2 Only in the non-interacting case at zero temperature all bosons can occupy the same 

single-particle state, </>(r), and therefore we may write the wave function of the N- particle 
system as ip( ri, r 2 ,.... rjv) = 0(ri)<?i>(r 2 ) ■ ... • <j>{ rjv) with normalization f d 3 rj</>(ri)| 2 = 1. 
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3 Split Operator method — a scalar version 

In this section, just for completeness, we briefly present SOM for a scalar order 
parameter (for more details see Refs. [2HI221[231[23] ) - Our goal is to find ip(t) 
satisfying m with the initial condition at t = 0 given by ip( 0). To this end 
we divide the time t into N intervals At such that t = N At. Then, assuming 
that At is small, we can write (with the first order of accuracy in At) 


N -1 


v>(t) = n exp (— iAtH(nAt))ip(0), 


(16) 


n —0 


where H(t ) is the right hand-side operator from ijCTl . In particular one has 


i/j(t + At) = exp 


1-72 i 1 2 


-i ( --V 2 + ^r 2 + g\i/j(t)\ 2 + V dip (t) ) At 


v>(t), (17) 


where for convenience only the dipolar contribution was denoted by V d i p . Using 
e A + B ~ e A e B from the Baker-Hausdorff theoren(f| we obtain 


i/j{t + At) ~ exp 


.At. 2 

*T V 


exp 


-iAt[-r + g\ip(t)\ + V dip (t) 


4>{t). (18) 


The above approximation is getting justified when At -C 1 (because commu¬ 
tator of the A and B operators appearing in the Baker-Hausdorff theorem is 
Zit-dependent, giving (At) 2 dependency in total). The right-hand side of Eq. 
(fT51) is just the Lie-Trotter splitting [55] applied to the GP equation (Hill . It is 
easy to improve the convergence of the above scheme by symmetrizing it with 
respect to one of the operators: A or B. It can be shown that the one-step 
formula of the form 


ip(t+At) ss exp 



exp 


-iAt. [ ^r 2 


g\i/j(t)\ 2 + V dip (t) 


exp 



(19) 

is of the second order in time. It is called the Strang splitting [55]. The operator 
exp(i(Z\f/4)V 2 ) does not depend on time. Therefore, when two successive time 
steps of the evolution scheme m are accomplished the neighbouring operators 
exp(i(Z\f/4)V 2 ) merge and take the form of exp(*(zlt/2)V 2 ) . 

Hence, effectively a single time-step At evolution (TT5T) can be split into two 
steps (only the first and the last steps in the sequence (fI51) should be treated 
differently). First, we define an auxiliary function 


V>i ( t) = exp 


-iAt I ^r 2 + g\%jj(t)\ 2 + V dip (t) 


^(t), 


( 20 ) 


3 The Baker-Hausdorff theorem reads e A+B = e A e B e 2 if [ A , B]] = 

0 and [B, [A, £?]] = 0. In our case the definitions of a A and B are A = i^-V 2 , B = 

-iAf (\r 2 + g\i>(t)\ 2 + V dip (t)) and are At- dependent, and e 2 oc I (because it’s 

cx At 2 ), so one can write e A + B ~ e A e B . 
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which we transform to the k-space using the Fourier transform 

$1 = (21) 

and then we evaluate the expression: 0 2 = exp[i^V 2 ]0i. In fc-space it be¬ 
comes a simple multiplication: 


■02 = exp 



( 22 ) 


Doing this, we compute the order parameter at time t + At, but it is still in 
momentum space. So, the last step is to return to the coordinate space (with 
the help of the inverse Fourier transform) 

0(t + At) = (23) 


4 Split Operator method — a spinor version 


In this section we introduce our extension to the SOM. For F = 1 atoms 
the order parameter 0 consists of three components 0 = (0_|_, 0o, 0-) which 
satisfies the equation similar to the Eq. m 


. d 


dt 



l -S7 2 + V(v,t)) ^0oj- 


(24) 


Here by V (r, t) we denoted the sum of the trapping potential and the nonlinear 
terms appearing in the right-hand side of Eq. (1151) describing the interactions 
between atoms. Then we use the scheme described in the previous section and 
obtain 


/ Vh 
V’o 
\ 0— 


(t + At) 



r.At 2 i 

r 

exp 


exp 


iAtV{r,t) 



(*)■ 


(25) 


Again, we introduce an auxiliary function 0i (this time as a 3-component 
vector), which we are going to transform to the momentum space. But first we 
need a careful treatment of a ‘potential term’ exp [— iAtV(r, t)] since V(r,t) 
is a matrix. The evolution in the position space requires calculation of the 
following expression: 

(M 

exp[-iAtV(r,t)] I 0° I (f). (26) 

To calculate (EUll we bring the V (r, t) matrix to the diagonal form: 


-l 


(X+ 0 0 \ 

o Ao 0 , 

V 0 0 A_/ 


1/ = PDP 


D = 


(27) 
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which, after utilizing the expression exp a: = J2 n gi yes 


exp[— iAtV] 


/V>+\ / 

I Vo I (t) = (l - iAtPDP- 1 


i (-iAtPDP~ 1 ) 2 + ..}j 



(28) 


The above expansion greatly simplifies after applying V n = PD n P 1 formula. 
One gets 


(t). 


exp [-iAtV] ^ ipo ^ (t) = P ^1 - iAtD + ^(- iAtD) 2 + ... ) P 1 ^ ip 0 j (*)■ 

(29) 

Next we can collect terms with D n back to the exponential form 

/ V>+ \ ( Vu \ 

exp [—iAtV] I ipo (t) = P exp[—iAtD]P 1 I V’o (t). (30) 

U-J U-) 


One can see now that we converted a problem of calculating (1M1) to calculation 
of the exponential function of a D matrix which is diagonal. So, we have 

/V>+\ / e ~ iX+At 0 \ /V>+\ 

exp [-iAtV] [ipo (i) = P 0 e ~ iX ° At o \ p -i ( t y 

\ip~ ) \ o o e - iX - At ) \V’-/ 

(31) 

The idea of calculating (l26l) by diagonalizing V(r,t) is the essence of our 
extension of the original SOkfl' A reader can think about (ED as a result of 
taking infinitely many terms in the Taylor expansion of the evolution operator 


/M . , 

exp [—iAtV] I V’o I (t) = |l — iAtV + 


{-iAtV) 2 + (-iAtV) 3 


2 ! 


3! 


( Va 
V’o 
\V- 


(t). 


(32) 

Utilizing (ED allows to overcome many problems: 1) we do not need to worry 
about the number of terms one should use to calculate right-hand side of (f32l) 
to achieve a good enough accuracy, and 2 ) diagonalization of small matrices 
(3 x 3 for F = 1 atoms or 7 x 7 for F = 3 atoms, like 52 Cr atoms) is significantly 
more efficient than calculating the right-hand side of (f3^1) , even assuming that 
only a few terms are taken into account. 


4 In fact we copied the original idea of the SOM: to calculate the evolution due to the 
kinetic part of the Hamiltonian one has to switch to the basis, in which the Laplace operator 
is diagonal, i.e., which is the momentum space. Because in the momentum space the dif¬ 
ferentiation becomes a simple multiplication by — zk, the evolution due to the kinetic term 
(which includes the second derivative in position space) becomes very easy to compute - 
instead of the second derivative one has to multiply Fourier components by exp[— k 2 At/2]. 
In the spinor version we have to do a similar trick in the position space: to facilitate calcu¬ 
lations of the evolution according to the potential energy, first we go to the basis in which 
the position dependent part of the Hamiltonian V(r,£) is diagonal, and only then we use 
the plane wave basis for evolution due to a kinetic part. 
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Once we successfully calculated the evolution due to the potential part 
of the Hamiltonian, we can follow remaining SOM steps as described in the 
previous section: first we use m to calculate an auxiliary spinor function 

/M 

0i = exp [~iAtV] I ip 0 I (t). (33) 

Next we move to the fc-space with -i/q 

■01 = -7 r [0i], (34) 


in order to calculate the evolution according to the kinetic energy part in the 
Hamiltonian 


02 = exp 



(35) 


Finally we come back to the coordinate space 


0(t + At) = T 1 [i>2\- 


(36) 


This sequence of steps propagates the initial wave function over the time in¬ 
terval At: ip(t) —> ip{t + At). 

It is clear from the above scheme, that SOM in a spinor version for F = 1 
atoms requires diagonalization of a 3 x 3 matrix at every spatial point on the 
grid and 3x2 = 6 Fourier transforms - a pair of forward and backward ones 
for each component of 0* (i goes from 1 to 3). 


5 Accuracy tests for spinor condensates 

We will present accuracy tests for real time propagation only, having already 
calculated the initial state by the imaginary-time propagation technique [26] 
applied to our extension to the SOM. We focus on real time propagation 
because the GP equation (fl5l) conserves the energy 

(E) = const , (37) 

the total norm 

N t ot = N + + N 0 + iV_ = const (38) 

with Ni = J|0i(r)| 2 d 3 r (i = +,0, —), and, assuming the dipolar interactions 
are neglected, the magnetization 


N + — N_ = const. (39) 

The expressions (l37l) , fl38[) and (l39l) are not any constraints for imaginary time- 
evolution. The energy decreases monotonically in every step of imaginary- 
time evolution till the ground state is reached [56]. The norm decreases as 
well and therefore after each single step of computations the wave function is 
normalized. Any algorithm attempting to solve the GP equation (fL5l) should 
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be verified with respect to how these constants of motion are preserved during 
the evolution. 

To test the conservation of the constants of motion (137H39I) (the case in¬ 
cluding the dipolar interactions will be discussed in the next section) we chose 
a system of 3 x 10 4 8 'Rb atoms, with a o = 5.387 nm and <22 = 5.313 nm (ac¬ 
cording to Ref. W) as contact interaction parameters, confined in a pancake 
trap (wj = ujy = 2 tt x 100 Hz, u> z = 2n x 2000 Hz). We start with all the 
atoms in rriF = 0 component and monitor the evolution of the system. Since 
the contact interactions m allow for the transfer of atoms from the initial 
state to the tof = ±1 states, one could expect the appearance of spin dy¬ 
namics. This dynamics was already investigated in [28] but here we will focus 
mainly on the numerical aspects, refereeing the reader to Ref. [28j for a deeper 
understanding of the physical background. 


iff 

~A 


LU 

V 



0.9998 


100 200 300 400 500 600 


1.0005 
1.0004 
^ 1.0003 
z 1-0002 
1.0001 
1.0000 

0 100 200 300 400 500 600 0 100 200 300 400 500 600 



TIME [r] 


TIME [r] 


TIME [t] 


Fig. 1 Illustration of the conservation of the total energy, the total norm, and the magne¬ 
tization of a spin-1 atoms system. Eq and Nq are the energy and the number of atoms at 
t = 0, respectively. Details of the system are given in the text. Colors depict different time 
steps used in simulations: At = 0.00025r (black), At = 0.0005r (blue), and At = O.OOlr 
(red). Time unit equals r = 1.59155 ms. 


We performed fully 3D simulations on a Cartesian grid of 2 5 2 5 2 4 points, 
with the spatial steps equal to Ax = Ay = 0.6, Az = 0.1 in oscillatory units 
(we have chosen Hco x as a unit of energy, thus the oscillatory unit of length 
equals a = 1.07961^im). Fig. [L] shows the conservation of the total energy, 
the total number of atoms, and the magnetization of a spin-1 gas for different 
time steps At. It is not surprising that the biggest time step leads to the worst 
conservation of the energy and the norm. This behavior is expected since At 
plays the crucial role in the SOM (compare At dependence of commutators 
in the Baker-Hausdorff theorem ©)• Even in the worst case presented here 
the energy is preserved with accuracy better than 0.1%, which is acceptable in 
many simulations. It is worth to notice that going into smaller At (for example 
decreasing it by a factor of 2) significantly improves accuracy (by a factor of 
2). It is important that we can easily control this accuracy by modifying At. 
Fig. [U also shows that the spinor SOM conserves the magnetization (l3iH) which 
is expressed here in absolute numbers. Let us remind that the total number 
of atoms in the system is N to t = 30000. 

Fig- [U shows a few interesting features: 1) conservation of the total energy 
starts to be worse at particular time (t ss 280 t), and 2) we note an almost 
linear behavior of ( E) /Eq curves - so one might be worried about conservation 
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of the energy for longer evolution time (t > 600r). To clarify the second 
issue we would like to remind the reader that we are dealing with ultracold 
atoms and that the typical lifetime of the condensate is of the order of a 
few seconds mmm- That is why there is no need to continue the evolution 
for more than here (1 s). Let us note however, that extending the total time 
by a factor of 5 (up to 5 seconds, which is a large evolution time from the 
point of view of real experiments) we get the total energy conserved with 1% 
accuracy (for At = 0.0005r) or even better - with 0.5% accuracy (in case of 
At = 0.00025r). That gives us an excellent level of conservation of the total 
energy, and moreover - we can control it by a proper choice of At. 

What happens around t r; 280t is shown in Fig. [2l left frame. Evidently, 
some nontrivial dynamics is triggered - atoms start to flow from tof = 0 to 
tof = ±1 components. Obviously, some spin dynamics might be expected in 
the system since our initial state is not the lowest energy solution of the spinor 
configuration. Spin changing collisions between atoms in the ttif = 0 state 
start to produce atoms in tuf = ±1 states. This reach spin dynamics evolves 
the system into the state of thermal equilibrium. This state is approached in 
a steplike process, see Ref. 1 281 . 


30000 

25000 

JI 20000 

J 15000 
10000 
5000 
0 


Ua> 

30000 
25000 
^ 20000 
^ 15000 


25000 

20000 
^ 15000 

~VWWWVWwv\ 


^ 10000 
5000 
0 

/Ay'W'vMn 

g 10000 

5000 

0 

JW/WWWWWV 


0 100 200 300 400 500 600 


0 100 200 300 400 500 600 


0 100 200 300 400 500 600 


TIME [r] 


TIME [r] 


TIME [r] 


Fig. 2 Time evolution of the number of atoms in each component for different values of 
initial seed. Left frame: N+( 0) = -/V_(0) = 10 -12 . Middle frame: iV_|_(0) = iV_(0) = 10 _1 . 
Right frame: 7V+(0) = 10 -12 , -/V_(0) = 4197. Other parameters: At = 0.0005r, 2 5 2 5 2 4 
points grid. Although a large transfer of atoms to mp = ±1 states starts at different times, 
the accuracy measured by a degree of conservation of the total energy and the total norm 
remains the same. For a deeper understanding of the physics behind we suggest to read Ref. 

m- 


Around time t ~ 280r a large number of atoms is transferred from initially 
populated ttif = 0 component to the other states (Fig. [2j left frame). One 
might wonder if this characteristic time at which the spin dynamics is triggered 
depends on the initial condition. As we have already mentioned, initially we 
assume N tot = 30000 atoms in tof = 0 component and almost zero atoms in 
ui f = ±1 states. To trigger the spin dynamics we need some seed in initially 
empty components, and, in fact, we used a seeding field in tof = ±1 states. 
The seed was implemented by choosing a complex random number at each 
point of the spatial grid. The seed plays a role of quantum fluctuations which 
are present in a real system. The quantum fluctuations are missing in the 
mean-field description (we have neglected all the quantum fluctuations, Eq. 
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(USD)- To see the spin dynamics @ the seed must be present at t = 0. We have 
checked, that changing the amount of seed (i.e. values of iV+(0) and AG(0) 
being the total number of atoms in mj? = +1, —1 states, respectively, at t = 0) 
results in shifting the time at which the significant transfer of atoms starts. For 
example for initial values IV+(0) = _/V_(0) = 10 -12 the spin dynamics begins 
around t ~ 280r (Fig. [2] left frame), but for N + ( 0) = 7V_(0) = 10 _1 it happens 
at t ps lOOr (Fig. [21 middle frame). We have also checked the imbalanced initial 
conditions by taking the initial populations as JV+(0) = 10 -9 and iV_(0) = 
4197. Fig. [21 (right frame) shows how the populations depend on time (brown 
solid line - N + (t), brown dotted line - The time at which the spin 

dynamics is triggered depends monotonically on the amount of the seed. What 
is important is that the value of the seed does not change qualitatively the 
dynamics - it changes only the time at which the non-trivial dynamics begins. 
Although the transfer of atoms starts at different times, the conservation of 
the total energy and the total norm (as well as the magnetization) is of the 
same order in all situations. 


6 Dipolar interactions 

It is well known that the dipolar interactions couple the spin and the orbital 
motion of colliding atoms. The projection of total spin of interacting atoms 
can change at most by 2 as implied by the expression (USD- When it happens 
it means that the spin goes to the orbital angular momentum of atoms. In 
other words, atoms changing their spin must acquire orbital motion. This is 
the famous Einstein-de Haas effect [20] which has been already discussed also 
for the systems of ultracold atoms It can be rigorously shown that 

the sum of the projections of the total spin and the total orbital angular 
momentum is preserved during collision. Indeed, the commutator [Vd,L ± z + 
L2z + Fiz + F2 Z ] equals zero, where L\ z and L2 Z denote the projections of 
the orbital angular momenta of colliding atoms [29] . Therefore, assuming the 
external potential has an axial symmetry along the z axis, the quantity L z +F z 
should be conserved during the evolution according to the GP equation (fl5l) 

L z + F z = const. (40) 

Any algorithm attempting to solve the GP equation (fl5l) with dipolar inter¬ 
actions should be verified against the quality of this constant of motion. 

But first we have to calculate the elements of Hd matrix (see the formula 
© and the discussion after it). One needs, actually, to obtain only two matrix 

5 By spin dynamics we understand the situation with non-negligible transfer of atoms from 
one component to the other. From Eq. (1151) it follows that '01 satisfies the equation i '((, = 

{Ho + H c Ii)-01 + H c ioV>o = {Ho + (co + C 2 )(|-01 | 2 |0ol 2 ) + (co - c 2 )|-0 —i| 2 )-0i + c 2 VA 1 i>oi>o 

and if 0i (0) = 0_i (0) = 0 the right hand side vanishes, i.e. the ipi field is constant what 

corresponds to a ‘spin frozen’ situation. To account for quantum fluctuations we need some 
small initial seeds in the fields 0i and 0_i. 
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elements, H d n 


H d n (r) =H 2 ^ 2 Jd 3 r' 


r — r 


./13 


- 3 


(*-*') 


M 21 


r — r 


'15 


x (IV’iO'Or — IV , -i( i " , )r) 


-3 

-3 


V2 

and H d io 


kl j rf3r ' [r - r'\5 [( x ~ x ') y ^\ X ^Wo( r ') + ^o( r ')ip- i( r ')) 

/ d3r ' |r — r '| 5 ^ ~ ^ + ^ ~ X «(r>i(r')+f-i(r>(r')) 

(41) 


ftV 


H dw (r) = -3^fd 


,[{x-x') - i(y - y')\{z - z') 


r'|5 


x (l^i(r')l 2 - |^-i(r')| 2 ) 


- 1 * 1 * f x W(r > 0 (r') + «(r>- l( r')) 


r — r 


+?i 2 7 2 / dV 


1 3 (a; - a;') 2 + {y - y') 2 


— r' 1 3 2 


r'|5 


X (V’o(r')V’i(r') + ^*i(r , )V'o(r')) • 

(42) 


All other nonvanishing elements can be expressed in terms of H d u and H d io- 

Integrals appearing in the expressions for the matrix elements (1411) and 
dH are the convolutions, therefore to calculate them it is convenient to use 
the Fourier transform technique. This is because the Fourier transform of the 
convolution is the product of the Fourier transforms of the functions which 
are convolved. The convolutions which contribute to (SB and HMD convolve 
functions which originate from the dipolar interactions with the ones com¬ 
posed based on the spinor components. Since the spinor wave function evolves 
according to the GP equation (SB, the Fourier transform of the part of SB 
and (1421) dependent on the spinor components is calculated numerically. On 
the other hand, the Fourier transform of the part which originates from the 
dipolar interaction is calculated analytically. The rest of this Section explains 
how it is done. 

We use the following convention of the Fourier transform m 


•^[/(r)] = /( P) = J d 3 r exp ^ + ^p ■ r j /( r). (43) 

To evaluate fl3l) it is convenient to transform coordinates in such a way that 
vectors p and z become parallel. This, of course, simplifies a scalar product 
p ■ r appearing in the argument of the exponential function. The required 
coordinate transformation is a product of two rotations: the first one is the 
rotation by an angle /3 around z axis and the second one by an angle a around 
already rotated y axis (see Fig. 0. The final rotation matrix is given by 


x 

y 

z 


( cos a cos /3 
cos a sin / 3 
— sin a 


— sin /3 sin a cos /3' 
cos /? sin a sin /? 

0 cos a 



cos a 


( 44 ) 
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Using this transformation we set the k vector (an argument of the Fourier 
transform, k = p /h) parallel to the z axis. Definitions of used angles, from 

Fig. [3] are: cosa = k z /k, sina = ^Jk 2 + k 2 /k, cos/3 = k x /y/k 2 — fc?, and 

sin/3 = ky/^/k 2 — k 2 . 



Fig. 3 The scheme of rotations used to calculate analytically Fourier transforms of functions 
which originate from dipolar interactions. 


kl 5 


As an example we will show in detail the calculation of the T jdp- — 3 
This particular Fourier transform gets important when the single component 
dipolar condensate, the system investigated already in the early days of BEC 
ESI, is considered. According to the definition C3l) we do integrate 


1 = 


/ 


d 3 r e 


ikr 



(45) 


with the help of the transformation of coordinates given by dm). Then, going 
to spherical coordinates and integrating over the azimuthal angle one obtains 

poo pi 

1= dr dt e lkrt — (l — 3 cos 2 a) (3t 2 — 1), (46) 

Jo J-i r 


where t = cos 8 and 8 is the polar angle. Now, doing double integration, first 
over t variable and then over the distance r we finally arrive at 

/ = — (l + 3cos(2a)) . (47) 

One has to be careful in calculating the expression (HH1) . It is easy to check 
that the result diverges if we first integrate over r variable. This is because, in 
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fact, the integral (H5l) does not converge and the regularization is required. It 
can be done based on the physical arguments saying that the size of the dipole 
is finite. Therefore, going back to (14611 one can first integrate it over r variable 
within the interval ( R , oo) getting 

I = J dt 7T (1 — 3 cos 2 a) (3f 2 — 1) T(0, — i k Rt ), (48) 


where T is the incomplete gamma function. Integrating over t variable and 
taking the limit R —> 0 lead us again to the result Gil). 

The other Fourier transforms we treat in a similar way and obtain m 


T 

T 

T 

T 

T 


4 - 3 ^ 


= —-ry-(l — 3 cos 2 a) 

O 


1 3 x 2 + y 2 


(x - iy)z 


l r l 

(a; + iy)z 


(x - iy) 2 


2tt o 

= — (1 — 3 cos 2 a) 
o 


= — 4 e 1/3 sin 2a 
o 

= — 4 e l/3 sin 2a 
o 

= 3 sin 2 a, 

O 


(49) 


From the numerical point of view calculations of Hdn and Hd io are the only 
change in m before we proceed to SOM algorithm. 


7 Accuracy tests for spinor dipolar condensates 

In Fig. H] we demonstrate how good is the SOM method with respect to the 
conservation of the sum of the total spin and the total orbital angular mo¬ 
mentum. We consider a system of N to t = 200000 87 Rb atoms confined in a 
spherically symmetric harmonic trap (ui = 2ir x 100 Hz). Initially all the atoms 
populate the ttjf = +1 Zeeman component. Then the magnetic field is turned 
on along the z axis. Provided the value of the magnetic field is resonant [29] , 
the atoms start, due to dipolar forces, to flow to other components (see Fig. 
a. As it is seen, significant number of atoms is transferred from the rriF = +1 
to ttif = 0 and mjr = -1 states. Certainly, the z projection of the total spin 
of the sample is changed during the evolution. However, Fig. [5] clearly shows 
that the projection of the total angular momentum is conserved very well. As 
in the original Einstein-de Haas effect the atoms in mj? =0 and mp = -1 
components start to rotate around the direction of the magnetic field showing 
vortices in tof = 0 and = — 1 components (see Fig. EH). 
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Fig. 4 Time evolution of the populations of all hyperfine states. Numerical parameters are 
as follows: Ax = Ay = Az = 0.5c*:, At = 0.0005r, grid size 2 5 2 5 2 5 , B z = —40 osc. units, 
fc 2 7 2 = 0.0000257722 osc. units. 



TIME [r] 


Fig. 5 Illustration of the conservation of the sum of the projections of the orbital angular 
momentum and the spin for various grid sizes. Figure shows the projection of the total 
angular momentum per atom in the system. Since we start calculations with all atoms 
being in the rrip = +1 Zeeman component, this quantity should be equal 1 all the time. 
Here we used At = 0.0005r and Ax = Ay = Az = Ar. 



0 


Fig. 6 Isodensity and phase of the order parameter at t = 0.14 s for the simulation presented 
in Fig. [4] Values of the density for mp = 1, 0, —1 (from left to right) equal 7.96 x 10 13 cm -3 , 
7.96 x 10 13 cm -3 , and 7.96 x 10 12 cm -3 . The color on the surface shows the phase of the 
order parameter for a given component (color scale on the right). Note characteristic vortex 
structures in the mp = 0 component (it is a single quantized vortex) and in the mp = — 1 
component (it is a doubly quantized vortex). 
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8 Degenerate Fermi gases 


Ultracold Fermi gases serve themselves as ideal quantum simulators because 
of highly controllable laboratory conditions they can be realized at. Two- 
component Fermi gases have been recently used to study the properties of 
strongly interacting fermionic systems SS- It is important to experimentally 
confirm the occurrence of a ferromagnetic instability driven by the short range 
repulsion, responsible for the transition from para- to ferromagnetic phase. 

It has been discovered already many years ago that the oscillations of 
the electron cloud in a many-electron atom can be viewed as a motion of a 
fluid characterized by the density and the velocity fields m- Such motion is 
described by the hydrodynamic equations [36]. Here we propose to use the 
hydrodynamic equations to study the properties of two-component fermionic 
gases m- Assuming the velocity fields are rotation-free, these equations read 

d 

—n± = -V(n± v±), 

d / 5T m 2 \ 

m Qi' r ± = -V ( + yV ± + Vtrap + 3 J , (50) 


where (nj( r, t), Vj(r, t)) denote the density and velocity fields of j —th (j = ±) 
component. T is the intrinsic kinetic energy of the gas and is calculated as 
in the Thomas-Fermi approximation [5] . Including the gradient corrections [8] 
one gets 


ST 2/3 , n 2 

5n± 1,1 2m y/n± 1 


(51) 


where £ = 1/9 and A = 6 5 / 3 ?i 2 7r 4 / 3 /(12m). Eqs. (I5U1) can be recast, by using 
the inverse Madelung transformation 1101 , to the pseudo-Schrodinger equation 


<4v ± = 


_^lv 2 + —(i_j)Z!l%l 
2 m 2rri |^±| 


^I4t| 4/3 + V trap 


1p± 


(52) 


Eqs. (l5^1l take the form of Eq. (liMl) . although here the system is described by 
two-component object (pseudo-wavefunction (//q_(r),^-(r)) 7 ) and therefore 
at each time step at each spatial point one has to diagonalize 2x2 matrix. 
Hence, the numerical method discussed in Sec. Q] can be directly used to study 
the dynamics of a Fermi mixture in the frame of Thomas-Fermi approximation. 

In Ref. [37] we investigated the ground state densities of repulsive two- 
component Fermi gases. Numerically, we just evolve the system according to 
Eqs. ([521) by using the imaginary time technique. By increasing the strength 
of repulsion we observe the transition from the identical density profiles for 
two species towards, first, isotropic and, finally, anisotropic separations of two 
components (see Fig. 5 in Ref. E2). This indicates indirectly the existence of 
a ferromagnetic instability in a system of two repulsive Fermi gases. 
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Here we investigate the accuracy of above discussed numerical algorithm 
after the trapping potential for two gases is periodically disturbed. In Fig. [7] we 
show how preserved is the total energy of the system after the trap is restored 
for different values of the repulsion strength. For the set of N + = N_ = 10 
atoms we analyze qualitatively different cases: the symmetric one correspond¬ 
ing to paramagnetic phase and the anisotropic separation case which repre¬ 
sents the ferromagnetic phase (see Fig. 5 in Ref. EH). The frequency of the 
trap is modulated for a short period as oj(1 + Asin(.!2i)) with the small (just 
to avoid any nonlinear effects) driving amplitude A = 5% and driving fre¬ 
quency fi = 2 uj (w is the unperturbed trap frequency). So, we investigate 
the monopole oscillating mode. After the trapping potential is restored we 
observe oscillations of both components. The insets in Fig. [7] show the oscilla¬ 
tions of the following quantities: f (x 2 + y 2 )n±(r)dr which are experimentally 
accessible after column density masurement along z direction is done. For fer¬ 
romagnetic phase the clouds oscillate with the frequency equals twice the trap 
frequency (the inset in the right frame in Fig. [7]). This is because in this phase 
the atoms of different spins practically do not interact with each other - the 
occupy opposite regions in the space (see Ref. EH). Therefore they behave as 
two independent ideal fermionic gases. On the other hand, for paramagnetic 
phase (left frame in Fig. [7]) the frequency increases (and for g = 7.0 equals 
2.08 w) because during the oscillations atoms of different spins interact all the 
time. 



TIME [i] TIME [i] 


Fig. 7 Illustration of the conservation of the total energy for a system of two-component 
Fermi gas. The number of atoms in each component is JV_|_ = N— = 10. The trap is initially 
disturbed (within the period up to the vertical dotted lines in each frame) by changing 
periodically the trap frequency. This pumps the energy into the system (main frames). 
When the trap is restored the energy is conserved at the level below one percent. The insets 
show the oscillations of each atomic cloud (both components behave in the same way) with 
the frequency which depends on the phase the system is in. The interaction strength changes 
from g = 7.0 (left frame) to g — 15.0 (right frame). 


Even more sophisticated description of a Fermi system which involves the 
single-particle spin-orbitals can be rewritten making possible to use the proce¬ 
dure detailed in Sec. [I] Let’s see how it works. We assume that the many-body 





























Unified way for computing dynamics of BEC and FG 


19 


wave function of N/2 + N/2 atoms is given by the single Slater determinant 


^(x 1, ...,Xjv) = 


Vm 


Pl(xi) ■ ■ 

. . ¥>i(xjv) 

<Av(x 1 ) . . 

. . ipn(x-n) 

for j = 1, 

.., N/2 and 


(53) 


V?j(x) = <Pj_ N/2 ( r )b{s) for j = N/2 + 1 which equally share the two 

spin states a(s) and b(s). The time-dependent Hartree-Fock equations for the 
spatial orbitals are then given by 


= U^ + 'Mr)) ^ 1} (r ,t) 

N/2 

+ J2 dr'\^\v/t)\ 2v aa {r-v') yf\r,t) 
k=i' 

N/2 

+ Y1 / dr ' l^fc 2) ( r, ^)| 2 Kb(r- r') ^(r.i) 
fc=i ' 

JV/2 . 

/ dr '(^i 1) ( r '^))* 1/ aa(r-r , )^ 1) (r',i) ^(M) (54) 

fc=i ' 

for j = l,...,iV/2, where the terms 14 a and 14b describe the interactions 
between atoms being both in the a(s) state and when one atom is in the a(s) 
state and the other in the b(s) one. Analogous set of equations is fulfilled by 
spatial orbitals yA (r,t), j = l,...,N/2. The first and the second integrals and 
the last one are called Coulomb and exchange terms, respectively. Equations 
(15411 for each kind of spatial orbitals can be written in the matrix form 

v {1) = (H 0 + yr + vs b - K“°) ^ (1) 

ih l. ^( 2 ) = ( Ho + V bb + yba _ ybb x) ^( 2 ) _ ( 55 ) 

where (</? (1) ) T = (^(r),<^} 2 (r)) T and (</? (2) ) T = (^ 2) (r),..., ^ 2 (r)) T . 
The Coulomb matrices are diagonal with equal elements, for example 
(VS a ) jk = J2if dr' lip^(r',t)\ 2 V aa {r - r ') S jk . The exchange matrices pos¬ 
sess off-diagonal elements, for example = f dr' (ip^ (r', t))* y, n (r — 

r') Lp'P (r', t). Surprisingly, the dynamics of a many-fermion system is again 
described by the equation like Eq. <ED in Sec. Q] Hence, the method intro¬ 
duced in Sec. [I] can be used. Here, however, at each spatial point we have to 
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diagonalize the square matrices of the size equal to the half of the number 
of atoms. Eqs. (GSD in the very simple case when only contact interactions 
between different spin states is allowed (the effective Hamiltonians on the 
right-hand side of (1551) 1 are then diagonal) was already studied by us in Ref. 
1381 . More demanding case when the pairing and ferromagnetic instabilities 
compete each other [23 (and the Hamiltonian matrices in (1551) are full) is 
under investigation 


9 Conclusions 

To summarize, we have presented the extended version of the SOM algorithm 
which turns out to be very efficient in simulating the evolution of spinor BEC 
systems with nonlocal interactions. We use this extension to solve the set of 
nonlinear partial integro-differential equations. In fact, this algorithm can be 
used to describe the evolution of any multicomponent system. It might be 
a spinor condensate of rubidium, chromium, erbium, or dysprosium atoms 
as well as the mixture of bosonic species. The algorithm can be also ap¬ 
plied to multicomponent systems consisting of indistinguishable or distinguish¬ 
able fermionic atoms. The combination of parallelization (with OpenMP tech¬ 
nique), the fast Fourier transform (via FFTW routine), and the linear algebra 
algorithm for diagonalization (e.g., from lapack packages) to the multicom¬ 
ponent SOM algorithm can be run very efficiently even on a typical desktop 
computer. We have proven that the algorithm conserves all constants of motion 
to a high accuracy. 
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